High pressure diamond-like liquid carbon 
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We report density-functional based molecular dynamics simulations, that show that, with increas- 
ing pressure, liquid carbon undergoes a gradual transformation from a liquid with local three- fold 
coordination to a 'diamond-like' liquid. We demonstrate that this unusual structural change is well 
reproduced by an empirical bond order potential with isotropic long range interactions, supple- 
mented by torsional terms. In contrast, state-of-the-art short-range bond-order potentials do not 
reproduce this diamond structure. This suggests that a correct description of long-range interactions 
is crucial for a unified description of the solid and liquid phases of carbon. 



When a simple fluid, such as argon, is cooled below the 
critical temperature, it can form a liquid phase. However, 
more complex fluids, such as molten phosphorus, can 
form two distinct liquid phases that are separated by a 
first-order liquid-liquid phase transition (LLPT)i. There 
is increasing evidence that a variety of network-forming 
liquids can undergo an LLPT (see e.g. Franzese et al^). 
In particular, it has long been suspected that elemental 
carbon may exhibit an LLPTiSi^. As direct experiments 
are difficult in the relevant temperature regime (4000 
to 6000 K), much of the evidence for such a transition 
comes from simulations. Classical simulations based on 
the Brenner bond-order potential with torsional terms^ 
predicted a coexistence line between a 2-fold and a four- 
fold coordinated liquid starting from the melting line of 
graphite end ending in a critical point (at ~ 8800 K)i 
as shown in Fig. ^ Recent first-principle investigations^ 
have ruled out the occurrence of such a transition in the 
range of specific volumes 6.59 - 15.6 A'^/atom at 6000 
K. More importantly, these ab-initio results found a ma- 
jority of three-fold coordination over the studied range. 
Wu et al.^ pointed out that the presence of 7r-bonds in 
the liquid phase makes it difficult to represent the tor- 
sion potential adequately. They concluded that empirical 
potentials are unable to capture these subtle electronic 
effects. 

In this Rapid Communication, we reexamine these is- 
sues by extending the range of studied densities while 
comparing density functional theory based molecular 
dynamics (DF-MD) simulations to Monte Carlo (MC) 
results based on two state-of-the-art empirical poten- 
tials. The first is the recently proposed long-range car- 
bon bond-order potential (LCBOP— ) extended with tor- 
sional interactions; the second is the reactive empirical 
bond-order (REBO^- ) potential which improves the ear- 
lier version due to Brenner^. LCBOP introduces long 
range terms to account for interplanar binding of graphite 
which is not described by the REBO potential. 

On the basis of DF-MD we find that, whilst no true 
LLPT occurs even for higher densities than in Wu et al&, 



a well defined switching of the dominant coordination 
from three- to four is found at ^ 5.8 A'^/atom, just be- 
yond the range studied in Wu et al.^. This 'diamond- like' 
four-fold coordinated liquid is highly structured, with 
a strongly anisotropic angular distribution of the first 
neighbors. Interestingly, the empirical potential LCBOP 
gives an impressive agreement with the DF-MD results 
for the structural properties of this high density liq- 
uid whereas the short-ranged REBO potentialiS yields 
a graphite-like liquid in this range. It should be noted 
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FIG. 1: (Color online) Equation of state of liquid carbon at 
T=6000K. The inset shows the P — T phase diagram, with di- 
amond (D), graphite (G) and liquid (L) domains: the dotted 
curve is the sampled isotherm; the solid curves are experi- 
mental data: the diamond-graphite coexistence curve from 
Bundy^i, the graphite melting curve from Togayai^. The 
dot-dashed line is the estimated diamond melting line from 
Glosli and Res—, obtained with Brenner potential; the dashed 
curves are the liquid-liquid coexistence curve ending in a criti- 
cal point and its graphite melting curve from Glosli and Rea-. 
The solid o in the inset indicates the switching in coordination 
shown in Fig. |21 
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that both the low-density three-fold coordinated and high 
density 'diamond-like' liquid are very unusual. In fact 
most covalent semiconductors become highly coordinated 
in the liquid phase Highly coordinated phases have 
been predicted by Grumbachi^ also for carbon albeit at 
pressures of TPa. 

We performed constant volume DF-MD simulations 
using the Car-Parrinello methodic as implemented in 
CPMD packageii. The system consisted of 128 atoms in 
a cubic box with periodic boundary conditions at 9 den- 
sities and a temperature T — 6000 K, imposed by means 
of a Nose-Hooveoi^ thermostat. We used the Becke^^ ex- 
change and Perdew^ii correlation gradient corrected func- 
tional (BP) with a plane wave basis set cut off at 35 Ry 
and sampled the Brillouin zone only in the gamma point. 
BP gives a correct description of bulk diamond. Each 
state point was studied for 5 ps, starting from a sample 
equilibrated via LCBOP; only minor structural changes 
occurred in the first tenths of ps. Since liquid carbon 
is metallic, we imposed a thermostat for the electronic 
degrees of freedom in order to ensure a proper imple- 
mentation of the Car-Parrinello scheme^^ . We performed 
the MC simulations of 128 particles in a cubic box with 
periodic boundary conditions with the LCBOP and the 
REBO potential. We sampled at 6000 K the constant vol- 
ume ensemble for specific volumes larger than 8 A^/atom 
and the constant pressure ensemble for smaller specific 
volumes where the increase of pressure is steeper, with 
an overlap between the two regions to check for consis- 
tency. MC simulations at selected volumes with 512 and 
4096 particles using the LCBOP showed negligible differ- 
ences with the 128 particle results for the local structure 
and equation of state. 

In order to reproduce the ab-initio results, we needed 
to refine the LCBOP. The detailed procedure will be de- 
scribed elsewhereSi. Here we give a brief summary of the 
changes. Following the strategy of the REBO potentialiS 
we have introduced a correction to the angular function 
in order to stabilize small clusters in shapes other than 
chain-like. It is in fact known (see e.g. Raghavachari 
and Binkley^'^) that, while odd numbered clusters prefer 
to arrange in a chain, even numbered ones can also have 
metastable ring or cage structures. We focused attention 
on a planar C4 cluster with C2V symmetry (a rhombus 
with two ~ 60° anglesS^) and a cubic Cg^. It turned out 
that the angular dependence of the bond order had to be 
weakened for coordination lower than three. More im- 
portant is the addition of torsional interactions, in such 
a way as to describe the conjugation dependence of the 
torsional energy barrier in agreement with the ab-initio 
calculations for double and conjugated bonds of Wu et 
al&. This dependence is not captured by the REBO po- 
tential. There, the barrier as a function of the dihedral 
angle Uijki is always described by sin'^{ujijki) and differs 
only by a scale factor between these two bonding sit- 
uations. In the original LCBOP, a variable N'^"^-' was 
already defined, assuming values 1 or respectively in 
these two extreme cases, to account for the presence or 
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FIG. 2; (Color online) Coordination fractions versus specific 
volume of liquid carbon. Note that aX V ~ 6 A'^/atom DF- 
MD and LCBOP predict a switch from three- to four-fold co- 
ordination whereas REBO evolves from a 2-fold to a graphite 
like liquid. 



absence of a tt bond. We fit the results of Wu et al& 
for N^j^-' = 0,1 by two polynomials in cos'^{uJijki) and 
interpolate between them by means of a function decay- 
ing rather quickly away from N^""'' — 1. The torsional 
energy is activated only for a bond between three-fold co- 
ordinated atoms, as for the REBO potential. The smaller 
torsional energy of a bond between four-fold coordinated 
atoms is already taken into account by the long range 
interactions between the third neighbours. 

In Fig. ^ we compare the pressure- volume isotherms 
of our DF-MD simulations with our MC results based 
on the empirical REBO potential and LCBOP, and with 
the DF-MD data of Wu et al.^, all at a temperature of 
6000 K. In view of the relatively low cut-off (35 Ry), we 
had to correct the pressures for the spurious contribution 
due to Pulay forcesS^. In the density range where we can 
compare with the results of Wu et al&, the pressures that 
we compute are some 15% lower than those reported by 
Wu et al. This relatively small difference is probably 
due to the different choice of density functionals. We 
have checked that our samples were indeed liquid. Over 
the whole isotherm we have observed diffusive behavior 
in both the MC-LCBOP and the DF-MD simulations, 
the latter indicating a self-diffusion coefficient at least of 
order lO^^cm^/s. 

Up to 60 GPa both empirical potentials seem to be in 
fair agreement with ab-initio data. However, the coordi- 
nation fractions shown in Fig. |21 indicate that the REBO 
potential fails to reproduce the correct liquid structure. 
It predicts, in the range between 6 and 7.5 A'^/atom, 
mainly 2-fold coordination whereas the present DF-MD 
results (and those of Wu et aim-), as well as the LCBOP, 
show dominant three-fold coordination and low 2-fold co- 
ordination. 

At higher pressures, DF-MD predicts a marked switch- 
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FIG. 3: (Color online) Pair correlation functions at 5000 K. 
Diamond (dotted curve) is at 150 GPa; the LCBOP liquid 
(solid curve) is at 150 GPa and 5.14 A^/atom; DF-MD liquid 
(dashed curve) is at 5.14 A^/atom. 

ing of dominant coordination from three to four around 
5.8 A^/atom. Judging from the isotherm of Fig. the 
transition seems to be continuous with no sign of a van 
der Waals loop. These results are consistent with the 
tight binding MD simulations of Morris et alm^. In con- 
trast, between 5.6 and 6.0 A^/atom, where the switch of 
dominant coordination takes place, the MC results based 
on LCBOP display large fluctuations in density at the 
imposed pressure of 100 GPa, resulting in a slight bend- 
ing of the isotherm of Fig^ Again, qualitatively the 
coordination fractions obtained by LCBOP agree with 
the DF-MD results even though with a more pronounced 
switch to coordination four. Curiously, the REBO po- 
tential predicts a first order phase transition in the same 
density range but to a completely different structure, a 
three-fold graphite-like liquid with well ordered sliding 
sheets which eventually get stuck upon further increas- 
ing of the pressure. The REBO potential never gives rise 
to four-fold coordination. 

If the three- to four-fold coordination change extrap- 
olates to lower temperatures, it might provide an expla- 
nation for the sudden change of the slope of the graphite 
melting linei^. The origin of the latter is still an open 
question, although it has been suggested that it is asso- 
ciated with a first order phase transition^ii^. Our results 
provide no evidence of a first-order transition but rather 
indicate a pronounced but continuous change of domi- 
nant coordination. 

Both the DF-MD and LCBOP simulations of the high 
density four-fold coordinated liquid show a strong dia- 
mond like positional and orientational order, as shown 
in Figs. EI and In Fig. |3| we compare the pair 
correlations g{r) of this high density liquid with that of 
a (meta)stable bulk diamond near the estimated coex- 
istence point, at 5000K and 150 GPa (also the LCBOP 
and DF-MD liquid samples were equilibrated at that tem- 
perature). One can see that up to the second neighbor 
shell the liquid has a structure almost as pronounced 




FIG. 4: (Color online) Angular correlation functions of first 
neighbors with state parameters as in Fig. |3 

as diamond. The g(r)'s obtained by LCBOP globally 
agree fairly well with the ones obtained by DF-MD—, 
except around 2 A, where the LCBOP minimum is too 
deep. In Fig. 0] we present the calculated angular correla- 
tion g^^") (9) for first neighbors, i.e. those atoms that fall 
within the short range cut-off of the LCBOP. Again, the 
first shell of neighbors in the liquid has a strong tetrahe- 
dral ordering, comparable to bulk diamond. To test how 
far the local diamond structure persists in the liquid, we 
define the angular correlation function for second neigh- 
bors (i.e. all those particles that are first neighbors of the 
first neighbors, excluding the central atom and the first- 
neighbor shell). Fig.[Slshows that, in the second-neighbor 
shell, the diamond structure is completely lost. Yet, the 
angular distribution is not structureless: we find a peak 
around 60° and a shoulder at ~ 35°; in a diamond lattice, 
the latter feature can be attributed to cross-correlations 
between the first and the second neighbor shells. 

It is rather surprising that the LCBOP potential re- 
produces the transformation to a predominantly four- 
fold coordinated liquid, while the REBO potential does 
not. Apparently, the isotropic long-ranged interactions 
play a crucial role. To demonstrate this, we verified that 
he short-ranged CBOP reproduces the behavior of the 
REBO potential^. This behavior is rather puzzling as 
long-ranged interactions are expected to play a negli- 
gible role at these high densities (see, e.g., Glosh and 
Reei^). Moreover, long-ranged interactions were intro- 
duced in Los and Fasolinc^ to describe three-fold coor- 
dinated graphitic phases, and no attempt was made to 
make the long range interactions dependent on the lo- 
cal environment. Torsional interactions appear to im- 
portant, since, without them, the calculated pressures 
would be too high for high densities and too low at low 
densities. 

We conjecture that the combination of torsional inter- 
actions, and long range forces is required to give the best 
description of the liquid. It would be interesting to check 
this conjecture by studying the liquid with other empir- 
ical potentials, in particular with Marks' environment- 
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FIG. 5; (Color online) Angular correlation functions for sec- 
ond neighbors (see text) with state parameters as in Fig. |21 
The peaks for the diamond are spread around the theoretical 
values for an fee lattice of -1, -0.5, 0, 0.5, with weights 6/66, 
24/66, 12/66, 24/66. 

dependent interaction potential (EDIP)^. The study of 
the liquid with this potential presented in MarksS^ is lim- 
ited to volumes of 6.2 A'^/atom where it gives coordina- 
tion fractions comparable to those obtained here by DF 
and by LCBOP and in Wu et al.^. 

In summary, on the basis of DF-MD simulations, we 
predict that liquid carbon should exhibit a continuous 



transformation from a three-fold to a mostly four-fold co- 
ordinated liquid at high pressure. This liquid has a strong 
local angular ordering, reminiscent of the diamond struc- 
ture. We have compared this finding with the results 
of MC simulations based on two empirical bond-order 
potentials with torsional terms, the short-range REBO 
potential and the long-range LCBOP. We find that the 
latter reproduces quite well the structure of both low- 
density and high-density liquid. The REBO simulations 
instead display a marked first order phase transition be- 
tween a mostly 2-fold and a mostly three-fold graphite 
fluid in this range. It is tempting to speculate that the 
existence of a four-fold coordinate liquid at high densi- 
ties will greatly facilitate the nucleation of diamonds from 
dense liquid carbon. 
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